##Opis danych i cel projektu Dane przedstawiają inforamcje o pacjentach z niewydolnością serca.Pacjenici zostali podzieleni na 4 kategorie: 1)Palacz, niewydolność śmiertelna 2)Palacz, niewydolność nie śmiertelna 3)Nie palacz, niewydolność śmiertelna 4)Nie palacz, niewydolność nie śmiertelna
Zmienne v1,..,v7 określają: 1)wiek 2)poziom kinazy kreatynowej we krwi (mcg/l) 3)procentowy poziom krwi opuszczającej serce przy każdym skurczu 4)poziom płytek krwi 5)pozim kreatyniny w surowicy (mg/dl) 6)poziom sodu w surowicy (mEq/l) 7)okres kontrolny (w dniach)
Źródło:https://doi.org/10.1186/s12911-020-1023-5 https://archive.ics.uci.edu/dataset/519/heart+failure+clinical+records
##Analiza danych jednowymiarowych
head(dane)
## Kategoria V1 V2 V3 V4 V5 V6 V7
## 1 3 75 582 20 265000 2 130 4
## 2 3 55 7861 38 263358 1 136 6
## 3 1 65 146 20 162000 1 129 7
## 4 3 50 111 20 210000 2 137 7
## 5 3 65 160 20 327000 3 116 8
## 6 1 90 47 40 204000 2 132 8
any(is.na(dane))
## [1] FALSE
Na początku zaprezentuję hiostogramy oraz podstawowe statystyki opisowe dla wszytkich zmiennych
"podsumowanie"
## [1] "podsumowanie"
apply(dane[,-1],2,summary)
## V1 V2 V3 V4 V5 V6 V7
## Min. 40.00000 23.0000 14.00000 25100 1.00000 113.0000 4.0000
## 1st Qu. 51.00000 116.5000 30.00000 212500 1.00000 134.0000 73.0000
## Median 60.00000 250.0000 38.00000 262000 1.00000 137.0000 115.0000
## Mean 60.83612 581.8395 38.08361 263358 1.41806 136.6254 130.2609
## 3rd Qu. 70.00000 582.0000 45.00000 303500 1.00000 140.0000 203.0000
## Max. 95.00000 7861.0000 80.00000 850000 9.00000 148.0000 285.0000
"odchylenie standardowe"
## [1] "odchylenie standardowe"
sort(apply(dane[,-1],2,sd), decreasing = TRUE)
## V4 V2 V7 V1 V3 V6
## 97804.236869 970.287881 77.614208 11.894809 11.834841 4.412477
## V5
## 1.043906
"wariancja"
## [1] "wariancja"
sort(apply(dane[,-1],2,var), decreasing = TRUE)
## V4 V2 V7 V1 V3 V6
## 9.565669e+09 9.414586e+05 6.023965e+03 1.414865e+02 1.400635e+02 1.946996e+01
## V5
## 1.089740e+00
"skośność"
## [1] "skośność"
sort(apply(dane[,-1],2,skewness), decreasing = TRUE)
## V2 V5 V4 V3 V1 V7 V6
## 4.4406886 4.2072326 1.4549746 0.5525927 0.4203739 0.1271606 -1.0428705
Największe wariancje mają zmienne v4 oraz v2, zmienne te mają też największe średnie oraz odcylenie standardowe. Wszystkie zmienne poza zmienną v6 mają dodatnią skośność.
for(i in 1:7) hist(dane[,i+1],main=colnames(dane)[i+1])
##Rozkład wielomodalny i obserwacje odstające
Rozkład normalny przypominają zmienne V6 oraz V4. Natomiast jeżeli chodzi o wielomodalnośc to przetestuje zmienne V3,V7,V2 oraz V5.
shapiro.test(dane$V6)
##
## Shapiro-Wilk normality test
##
## data: dane$V6
## W = 0.93903, p-value = 9.215e-10
shapiro.test(dane$V4)
##
## Shapiro-Wilk normality test
##
## data: dane$V4
## W = 0.91151, p-value = 2.883e-12
Żadna zmienna nie ma rokładu normalnego
silverman.test(dane$V3,k=1)
## Silvermantest: Testing the null hypothesis that the number of modes is <= 1
## The resulting p-value is 0.04004004
silverman.test(dane$V7,k=1)
## Silvermantest: Testing the null hypothesis that the number of modes is <= 1
## The resulting p-value is 0.001001001
silverman.test(dane$V2,k=1)
## Silvermantest: Testing the null hypothesis that the number of modes is <= 1
## The resulting p-value is 0.07007007
silverman.test(dane$V5,k=1)
## Silvermantest: Testing the null hypothesis that the number of modes is <= 1
## The resulting p-value is 0.04204204
p-value poniżej 5% mają zmienne V3,V7,V5. Oznacza to, że mają roakłady wielomodalne. Za pomocą silvermanplot sprawdzę ilo modalne są te zmienne
silverman.plot(dane$V3)
silverman.plot(dane$V7)
silverman.plot(dane$V5)
Zwykresu zmienne V3 oraz V5 są trójmodalne natomiast zmienna V7 jest cztero lub pięcio modalna.
Teraz korzystając z testu Grubbsa sprawdzimy obserwacje odstające.
apply(dane[,-1],2,function(x) grubbs.test(x))
## $V1
##
## Grubbs test for one outlier
##
## data: x
## G = 2.87217, U = 0.97222, p-value = 0.576
## alternative hypothesis: highest value 95 is an outlier
##
##
## $V2
##
## Grubbs test for one outlier
##
## data: x
## G = 7.5021, U = 0.8105, p-value = 4.647e-13
## alternative hypothesis: highest value 7861 is an outlier
##
##
## $V3
##
## Grubbs test for one outlier
##
## data: x
## G = 3.54178, U = 0.95776, p-value = 0.05195
## alternative hypothesis: highest value 80 is an outlier
##
##
## $V4
##
## Grubbs test for one outlier
##
## data: x
## G = 5.99812, U = 0.87887, p-value = 9.134e-08
## alternative hypothesis: highest value 850000 is an outlier
##
##
## $V5
##
## Grubbs test for one outlier
##
## data: x
## G = 7.26305, U = 0.82239, p-value = 3.95e-12
## alternative hypothesis: highest value 9 is an outlier
##
##
## $V6
##
## Grubbs test for one outlier
##
## data: x
## G = 5.35423, U = 0.90348, p-value = 6.145e-06
## alternative hypothesis: lowest value 113 is an outlier
##
##
## $V7
##
## Grubbs test for one outlier
##
## data: x
## G = 1.99370, U = 0.98662, p-value = 1
## alternative hypothesis: highest value 285 is an outlier
Zmienne V1,V3,V7 nie mają obserwacji odstających
##Macierz korelacj, redukcja wymiaru
Macierz Pearsona
heatmaply_cor(cor(dane[,-1]))
Ogólnie dane są ze sobą słabo skorelowane. W delikatnym stopniu można zobaczyć korelację danych V1 i V5, oraz V6iV3.
heatmaply_cor(cor(dane[,-1],method = 'kendall'))
Macierz korelacji Kendalla jest bardzo podobna do macierzy Pearsona. Można na niej zobaczyć jednak odrobinę słabsze skorelowanie danych poza danymiV1 i v5 które są ze sobą skorelowane minimalnie bardziej.
prcomp(dane[,-1])->pca.dane
summary(pca.dane)
## Importance of components:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7
## Standard deviation 9.780e+04 9.7e+02 77.66 12.1 11.23 4.305 1.009
## Proportion of Variance 9.999e-01 1.0e-04 0.00 0.0 0.00 0.000 0.000
## Cumulative Proportion 9.999e-01 1.0e+00 1.00 1.0 1.00 1.000 1.000
df.pca<-data.frame(pc1=pca.dane$x[,1],pc2=pca.dane$x[,2],pc3=pca.dane$x[,3],kl=as.factor(dane$Kategoria))
plot(x=df.pca[,1],y=df.pca[,2], col=as.factor(dane[,1]))
plot_ly(df.pca,x=~pc1,y=~pc2,z=~pc3,color=~kl,type='scatter3d')
## No scatter3d mode specifed:
## Setting the mode to markers
## Read more about this attribute -> https://plotly.com/r/reference/#scatter-mode
sort(pca.dane$rotation[,1])
## V1 V5 V6 V7 V3
## -6.359893e-06 -4.606192e-07 2.802784e-06 8.343479e-06 8.733849e-06
## V2 V4
## 2.427180e-04 1.000000e+00
sort(pca.dane$rotation[,2])
## V2 V6 V5 V4 V3
## -9.999990e-01 -2.640333e-04 2.199171e-05 2.427134e-04 5.596420e-04
## V7 V1
## 7.733495e-04 9.861072e-04
sort(pca.dane$rotation[,3])
## V7 V3 V6 V2 V4
## -9.993492e-01 -6.271944e-03 -5.006069e-03 -7.403701e-04 8.810713e-06
## V5 V1
## 1.866727e-03 3.510970e-02
Pierwsza składniowa zawiera prawie 100% wariancji.PC2 jest zdominowane przez zmienną V2, PC3 przez zmienną V7.
##T-sne
dane.tsne<-Rtsne(dane[-1],theta=0,dims=3,perplexity=5)
dane.tsne.2<-Rtsne(dane[,-1],theta=0.5,dims=3,perplexity=20)
summary(dane.tsne$Y)
## V1 V2 V3
## Min. :-58.489 Min. :-48.003 Min. :-48.013
## 1st Qu.: -8.766 1st Qu.:-22.906 1st Qu.:-17.219
## Median : -1.094 Median : 3.105 Median : -1.518
## Mean : 0.000 Mean : 0.000 Mean : 0.000
## 3rd Qu.: 15.767 3rd Qu.: 12.119 3rd Qu.: 17.640
## Max. : 64.226 Max. : 53.072 Max. : 54.333
summary(dane.tsne.2$Y)
## V1 V2 V3
## Min. :-11.9476 Min. :-18.7120 Min. :-13.299
## 1st Qu.: -6.8698 1st Qu.:-11.6139 1st Qu.: -9.430
## Median : 0.1425 Median : -0.5326 Median : -2.878
## Mean : 0.0000 Mean : 0.0000 Mean : 0.000
## 3rd Qu.: 6.9093 3rd Qu.: 10.4147 3rd Qu.: 9.972
## Max. : 10.4579 Max. : 21.4740 Max. : 19.852
dane.tsne.df<-as.data.frame(dane.tsne$Y)
dane.tsne.2.df<-as.data.frame(dane.tsne.2$Y)
plot_ly(dane.tsne.df, x = ~V1, y = ~V2, z = ~V3,color = dane$Kategoria, type = 'scatter3d')
## No scatter3d mode specifed:
## Setting the mode to markers
## Read more about this attribute -> https://plotly.com/r/reference/#scatter-mode
plot_ly(dane.tsne.2.df, x = ~V1, y = ~V2, z = ~V3,color = dane$Kategoria, type = 'scatter3d')
## No scatter3d mode specifed:
## Setting the mode to markers
## Read more about this attribute -> https://plotly.com/r/reference/#scatter-mode
Dane są ułożone chatycznie. Ze wzrostem parametru perplexity zbijają się w grupy, a poszególne grupy oddalają się od siebie.
##Analiza skupień Zacznę od ilości klasrów równej ilości kategorii
kmeans(dane[,-1],centers=4)->km.dane.2
df.pca<-data.frame(pc1=pca.dane$x[,1],pc2=pca.dane$x[,2],pc3=pca.dane$x[,3],kl=as.factor(km.dane.2$cluster))
plot_ly(df.pca,x=~pc1,y=~pc2,z=~pc3,color=~kl,type='scatter3d')
## No scatter3d mode specifed:
## Setting the mode to markers
## Read more about this attribute -> https://plotly.com/r/reference/#scatter-mode
Wygląda na to, że k-średnich ładnie podzieliło wykres, ale upewnimy się jeszcze
table(dane$Kategoria, km.dane.2$cluster)
##
## 1 2 3 4
## 1 5 6 3 16
## 2 17 16 3 30
## 3 19 17 2 28
## 4 20 31 6 80
mogło by być lepiej
wss <- rep(NA,9)
for(i in 1:9) wss[i] <- kmeans(dane[,-1], centers=i+1)$tot.withinss
plot(1:9,wss)
Z wykresu wynika, że 3 to optymalna liczba klastrów. Przeprowadzę więc klastrowanie dla tej liczby
kmeans(dane[,-1],centers=3)->km.dane.3
df.pca.3<-data.frame(pc1=pca.dane$x[,1],pc2=pca.dane$x[,2],pc3=pca.dane$x[,3],kl=as.factor(km.dane.3$cluster))
plot_ly(df.pca.3,x=~pc1,y=~pc2,z=~pc3,color=~kl,type='scatter3d',mode='markers')
Znowu wygląda ładnie, ale sprawdzimy to
table(dane$Kategoria, km.dane.3$cluster)
##
## 1 2 3
## 1 10 16 4
## 2 27 32 7
## 3 30 29 7
## 4 44 80 13
Znowu nie wygląda to dobrze, najwięcej wyników wrzuciło do klastra 3
Przeprowadze teraz klastrowanie aglomeracyjne oraz podziałowe.
plot(agnes(dist(dane[,-1], method='minkowski',p=1),diss=T))
plot(diana(dane[,-1]))
Wyniki są zbliżone ##Klasyfikacja pod nadzorem- lasy losowe
rfcv(trainx = dane[,-1],trainy = as.factor(dane$Kategoria))->rf.dane
rf.dane
## $n.var
## [1] 7 4 1
##
## $error.cv
## 7 4 1
## 0.4916388 0.5050167 0.5451505
##
## $predicted
## $predicted$`7`
## [1] 3 3 3 3 3 3 3 3 3 1 3 3 3 3 4 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 1 1 3 3 3 3
## [38] 3 3 3 3 3 3 3 1 3 1 3 3 4 1 1 3 4 3 3 3 3 3 1 3 3 3 3 3 3 1 3 3 3 3 2 2 2
## [75] 1 3 4 4 4 2 4 4 3 3 4 4 4 4 4 2 4 4 4 3 4 4 4 2 4 4 4 4 4 4 4 4 4 2 4 4 4
## [112] 4 3 4 4 4 4 3 4 3 3 4 4 4 3 4 3 4 2 4 4 3 4 4 4 4 4 3 4 2 4 2 4 4 4 4 4 4
## [149] 3 4 4 4 4 4 4 3 2 4 3 4 4 4 4 2 4 4 4 3 4 4 4 2 2 2 4 4 4 4 2 4 4 1 4 3 4
## [186] 3 4 1 4 4 3 4 4 4 3 4 4 2 4 4 4 4 4 3 4 4 4 3 4 4 4 4 4 4 4 4 4 2 4 2 4 4
## [223] 4 4 4 4 4 2 3 4 4 4 4 4 4 2 4 2 4 4 4 4 4 4 2 4 4 4 4 4 3 4 4 4 2 4 4 4 4
## [260] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 3 3 4 4 4 4 4 4 4 4 4 4 4 4 2
## [297] 2 4 4
## Levels: 1 2 3 4
##
## $predicted$`4`
## [1] 1 1 3 3 3 3 1 3 1 1 3 3 4 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 1 3 3 3 3 3
## [38] 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 4 3 3 2 3 3 1 3 3 3 2 2 3 1 3 1 3 2 2 2 2
## [75] 1 3 4 4 2 2 2 4 3 3 4 4 4 4 4 3 4 4 4 4 4 4 4 3 4 4 4 4 4 4 4 3 4 2 4 4 4
## [112] 4 3 4 4 4 4 3 4 4 4 4 4 2 4 4 3 4 4 2 4 3 4 4 3 4 4 3 4 2 3 2 4 4 3 4 4 4
## [149] 3 4 4 2 4 4 4 1 4 1 4 4 4 4 4 2 2 3 4 3 2 4 4 2 4 3 4 4 4 3 2 4 2 1 4 3 1
## [186] 3 4 1 4 4 3 2 4 4 3 4 4 2 4 4 4 4 4 3 4 4 4 4 4 4 4 4 4 3 4 4 4 2 4 4 3 4
## [223] 4 4 3 4 4 4 3 2 4 4 2 4 4 2 2 2 4 4 2 2 2 4 4 4 4 4 4 4 2 2 4 4 2 4 4 4 4
## [260] 4 4 4 4 4 4 4 3 4 4 4 4 4 4 4 4 4 4 4 2 4 4 4 2 4 4 4 2 4 4 4 4 2 4 4 2 2
## [297] 2 2 2
## Levels: 1 2 3 4
##
## $predicted$`1`
## [1] 1 1 3 1 1 3 1 1 1 1 1 1 3 3 3 3 3 3 3 3 3 1 3 1 4 4 1 3 1 3 1 3 3 3 4 4 4
## [38] 4 3 3 3 1 3 3 1 1 1 3 1 3 3 3 3 3 3 2 2 2 2 3 3 2 3 3 1 4 1 3 3 3 3 1 2 1
## [75] 1 1 4 2 2 2 2 4 2 3 4 4 4 4 4 4 4 2 3 4 4 4 4 4 4 4 4 2 4 4 2 2 1 2 1 4 2
## [112] 3 1 2 2 2 4 4 4 4 4 4 4 4 4 3 4 2 4 2 2 2 4 4 4 2 4 4 4 4 4 4 4 4 4 2 2 3
## [149] 4 3 4 2 4 2 4 4 4 4 4 4 2 4 2 2 3 3 3 2 3 2 4 2 2 2 2 2 4 2 2 4 4 4 1 1 3
## [186] 3 3 3 3 4 4 4 4 3 2 2 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 2 2 2 3 4 3 4 4 4 4 2
## [223] 4 2 4 2 4 2 3 4 4 4 2 4 4 2 4 4 4 4 2 4 2 4 4 4 4 4 4 4 4 4 4 4 4 4 2 2 2
## [260] 4 4 4 4 4 4 1 2 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 2 2 2 2 4 4 4 2 4 2 2 2
## [297] 2 2 2
## Levels: 1 2 3 4
Błąd walidacji krzyżowej jest najniższy przy 4 cechach
randomForest(dane[,-1],as.factor(dane$Kategoria),importance=T)->rf.dane.imp
varImpPlot(rf.dane.imp)
4 najlepsze cechy to V7,V5,V3,V1
randomForest(dane[,c(8,6,4,2)],as.factor(dane$Kategoria))->rf.dane.sel
rf.dane.sel$confusion
## 1 2 3 4 class.error
## 1 7 1 19 3 0.7666667
## 2 1 11 5 49 0.8333333
## 3 10 5 37 14 0.4393939
## 4 1 28 17 91 0.3357664
1-mean(rf.dane.sel$confusion[,5])
## [1] 0.4062099
Używając tej metody osiągamy średnią dokładność równą 42%. Największy błąd jest w klasie 2 i wynosi 82%
##KNN Zacznę od stworzenia danych testowych i treningowych
podział<-createDataPartition(dane[,1],p=0.8,list=F)
dane.train<-(dane[podział,])
dane.test<-(dane[-podział,])
prop.table(table(dane$Kategoria))
##
## 1 2 3 4
## 0.1003344 0.2207358 0.2207358 0.4581940
prop.table(table(dane.train$Kategoria))
##
## 1 2 3 4
## 0.1083333 0.2125000 0.2208333 0.4583333
prop.table(table(dane.test$Kategoria))
##
## 1 2 3 4
## 0.06779661 0.25423729 0.22033898 0.45762712
Teraz zajmę się skalowaniem i standaryzacją
sis_tr<-preProcess(dane.train[,-1])
sis_tr
## Created from 240 samples and 7 variables
##
## Pre-processing:
## - centered (7)
## - ignored (0)
## - scaled (7)
Żadna cecha nie została pomienięta
dane.train.sc<-predict(sis_tr,dane.train)
apply(dane.train.sc[,-1],2,sd)
## V1 V2 V3 V4 V5 V6 V7
## 1 1 1 1 1 1 1
summary(dane.train.sc[,-1])
## V1 V2 V3 V4
## Min. :-1.78861 Min. :-0.595577 Min. :-2.014368 Min. :-2.50869
## 1st Qu.:-0.77539 1st Qu.:-0.488644 1st Qu.:-0.672388 1st Qu.:-0.53879
## Median :-0.09991 Median :-0.344352 Median :-0.001398 Median :-0.02151
## Mean : 0.00000 Mean : 0.000000 Mean : 0.000000 Mean : 0.00000
## 3rd Qu.: 0.74443 3rd Qu.: 0.009746 3rd Qu.: 0.585718 3rd Qu.: 0.42149
## Max. : 2.85530 Max. : 7.719756 Max. : 2.682562 Max. : 5.09820
## V5 V6 V7
## Min. :-0.3936 Min. :-5.16093 Min. :-1.6628
## 1st Qu.:-0.3936 1st Qu.:-0.57020 1st Qu.:-0.7442
## Median :-0.3936 Median : 0.08562 Median :-0.1566
## Mean : 0.0000 Mean : 0.00000 Mean : 0.0000
## 3rd Qu.:-0.3936 3rd Qu.: 0.74144 3rd Qu.: 0.9471
## Max. : 7.0154 Max. : 2.49029 Max. : 1.8949
Wartość średnia i odchylenie standardowe są równe 0 i 1
set.seed(length(dane.train))
dane.test.sc<-predict(sis_tr,dane.train)
#jako k przjmę pierwiastek z liczby obserwacji
k<-sqrt(length(dane.test.sc$Kategoria))
dane.knn<-knn(dane.train.sc[,-1],dane.test.sc[,-1],dane.train.sc$Kategoria,k)
table(dane.knn,dane.test.sc$Kategoria)
##
## dane.knn 1 2 3 4
## 1 5 0 3 2
## 2 4 5 1 1
## 3 6 3 20 4
## 4 11 43 29 103
mean(dane.knn==dane.test.sc$Kategoria)
## [1] 0.5541667
Knn nie poradziło sobie najlepiej, średnia dokładność jest róna około 50%.
##Wnioski końcowe (1)Zrobiłem analizę podstawowych statystyk opisowych takich jak wariancja, skośność czy wielomodalność (2)Zanalizowałem macierz korelacji Pearsona, zastosowałem PCA i T-sne do redukcji wymiaru (3)Przeprowadziłęm analizę skupień za pomocą metdy k-średnich, aglomeracyjnej oraz podziałowej (4)Przeprowadziłem kalsyfikację z użyciem lasów losowych (5)Przeprowadziłem klasyfikację za pomocą algorytmu k- najbliższych sąsiadów
Pomimo zastosowania wielu metod statystycznych, analizy danych i
klasyfikacji stwierzdam, że wyniki króre otrzymałem pokazują, że
modelowanie tych danych może być trudne i mało skuteczne. Żeby osiągnąć
zadowalające efekty potrzeba byłoby więcej danych. Trzeba również
pamiętać, że używałem danych dotyczących pacjentów z niewydolnością
serca, a dane medyczne często nie są możliwe do
jednoznacznej interpretacji.